pacman::p_load(sf, terra, spatstat,
tmap, rvest, tidyverse)Hands-on -Exercise 4
Spatial Point Pattern Analysis (SPPA) is the evaluation of the pattern or distribution of a set of points on a surface. The points may represent:
events such as crimes, traffic accidents, or disease onsets, or business services (e.g., coffee shops and fast-food outlets) or facilities such as childcare centres and eldercare centres. First-order Spatial Point Pattern Analysis (1st-SPPA) focuses on understanding the intensity or density of points across a study area. It examines how the distribution of points varies over space, essentially identifying trends or patterns in point density. This type of analysis deals with the individual locations of points and their distribution, without considering interactions between them.
In essence, 1st-SPPA helps answer questions such as:
Where are points most densely located within the study area? Is point density uniform, or does it vary across space? How spread out is the point pattern? In this chapter, you will gain hands-on experience with spatstat to perform two commonly used 1st-SPPA methods:
The specific questions we would like to answer are as follows:
Are the childcare centres in Singapore randomly distributed throughout the country? If the answer is not, then the next logical question is where are the locations with higher concentration of childcare centres?
Installing and Loading the R packages
Importing and Wrangling Geospatial Data Sets
mpsz_sf <- sf::st_read("../Hands_on_Ex02/Data/geospatial/MasterPlan2019SubzoneBoundaryNoSeaKML.kml") |>
sf::st_zm(drop = TRUE, what = "ZM") |>
sf::st_transform(crs = 3414)Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source
`C:\Users\kchan\Desktop\chandru-ko\ISSS608-VAA\Hands-on_Ex\Hands_on_Ex02\Data\geospatial\MasterPlan2019SubzoneBoundaryNoSeaKML.kml'
using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84
Next, build a function called extract_kml_field for extracting REGION_N, PLN_AREA_N, SUBZONE_N, SUBZONE_C from Description field by using the code chunk below.
extract_kml_field <- function(html_text, field_name) {
if (is.na(html_text) || html_text == "") return(NA_character_)
page <- read_html(html_text)
rows <- page %>% html_elements("tr")
value <- rows %>%
keep(~ html_text2(html_element(.x, "th")) == field_name) %>%
html_element("td") %>%
html_text2()
if (length(value) == 0) NA_character_ else value
}mpsz_sf <- mpsz_sf %>%
mutate(
REGION_N = map_chr(Description, extract_kml_field, "REGION_N"),
PLN_AREA_N = map_chr(Description, extract_kml_field, "PLN_AREA_N"),
SUBZONE_N = map_chr(Description, extract_kml_field, "SUBZONE_N"),
SUBZONE_C = map_chr(Description, extract_kml_field, "SUBZONE_C")
) %>%
select(-Name, -Description) %>%
relocate(geometry, .after = last_col())mpsz_cl <- mpsz_sf %>%
filter(SUBZONE_N != "SOUTHERN GROUP",
PLN_AREA_N != "WESTERN ISLANDS",
PLN_AREA_N != "NORTH-EASTERN ISLANDS")write_rds(mpsz_cl,
"Data/mpsz_cl.rds")Next, code chunk below will be used to import the childcare Service data downloaded from data.gov.sg into R environment as sf data frame called chilcare_sf. The Simple Feature Geometry (sfg) of this geospatial data layer if
childcare_sf <- st_read("Data/ChildCareServices.kml") %>%
st_zm(drop = TRUE, what = "ZM") %>%
st_transform(crs = 3414)Reading layer `CHILDCARE' from data source
`C:\Users\kchan\Desktop\chandru-ko\ISSS608-VAA\Hands-on_Ex\Hands_on_Ex04\Data\ChildCareServices.kml'
using driver `KML'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84
Mapping the geospatial data sets
After checking the referencing system of each geospatial data data frame, it is also useful for us to plot a map to show their spatial patterns.
tmap_mode('view')ℹ tmap mode set to "view".
tm_shape(childcare_sf)+
tm_dots()Registered S3 method overwritten by 'jsonify':
method from
print.json jsonlite
tmap_mode('plot')ℹ tmap mode set to "plot".
Geospatial Data wrangling
Converting sf data frames to ppp class
spatstat requires the point event data in ppp object form. The code chunk below uses [as.ppp()] of spatstat package to convert childcare_sf to ppp format.
childcare_ppp <- as.ppp(childcare_sf)Next, class() of Base R will be used to verify the object class of childcare_ppp.
class(childcare_ppp)[1] "ppp"
You can take a quick look at the summary statistics of the newly converted ppp object by using the code chunk below.
summary(childcare_ppp)Marked planar point pattern: 1925 points
Average intensity 2.417323e-06 points per square unit
Coordinates are given to 11 decimal places
Mark variables: Name, Description
Summary:
Name Description
Length:1925 Length:1925
Class :character Class :character
Mode :character Mode :character
Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
(33590 x 23700 units)
Window area = 796335000 square units
Creating owin object
When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.
The code chunk below, as.owin() of spatstat is used to covert mpsz_sf into owin object of spatstat.
sg_owin <- as.owin(mpsz_cl)Again, class() will be used to verify the object class of sg_owin.
class(sg_owin)[1] "owin"
The result above confirmed that sg_owin is indeed in owin object.
sg_owin object can be displayed by using plot() function.
plot(sg_owin)
Combining point events object and owin object
In this last step of geospatial data wrangling, we will extract childcare events that are located within Singapore by using the code chunk below.
childcareSG_ppp = childcare_ppp[sg_owin]Nearest Neighbor Analysis (NNA) is a spatial statistics method that calculates the average distance between each point and its closest neighbor to determine if a pattern of points is clustered, dispersed, or randomly distributed.
Clark-Evans test is a specific statistical method used within NNA to quantify whether a point pattern is clustered, random, or uniformly spaced, using the Clark-Evans aggregation index (R) to describe this pattern. NNA provides a numerical value that describes the degree of clustering or regularity, and the Clark-Evans test calculates a specific index (R) for this purpose
Perform the Clark-Evans test without CSR
clarkevans.test() of spatstat.explore package support two Clark-Evans test, namely: without CRS and with CRS. In the code chunk below, Clark-Evans test without CSR method is used.
clarkevans.test(childcareSG_ppp,
correction="none",
clipregion="sg_owin",
alternative=c("clustered"))
Clark-Evans test
No edge correction
Z-test
data: childcareSG_ppp
R = 0.53532, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
Perform the Clark-Evans test with CSR
In the code chunk below, the argument method = “MonteCarlo” is used. In this case, the p-value for the test is computed by comparing the observed value of R to the results obtained from nsim (i.e. 39, 99, 999) simulated realisations of Complete Spatial Randomness conditional on the observed number of points.
clarkevans.test(childcareSG_ppp,
correction="none",
clipregion="sg_owin",
alternative=c("clustered"),
method="MonteCarlo",
nsim=99)
Clark-Evans test
No edge correction
Monte Carlo test based on 99 simulations of CSR with fixed n
data: childcareSG_ppp
R = 0.53532, p-value = 0.01
alternative hypothesis: clustered (R < 1)
Kernel Density Estimation Method
Kernel Density Estimation (KDE) is a valuable tool for visualising and analyzing first-order spatial point patterns. It is widely considered a method within Exploratory Spatial Data Analysis (ESDA) because it’s used to visualize and understand spatial data patterns by transforms discrete point data (like locations of childcare service, crime incidents or disease cases) into continuous density surfaces that reveal clusters and variations in event occurrences, without making prior assumptions about data distribution. It helps to begin understanding data distribution, identify hotspots, and explore relationships between spatial variables before performing more rigorous analysis.
Working with automatic bandwidth selection method The code chunk below computes a kernel density by using the following configurations of density() of spatstat:
bw.diggle() automatic bandwidth selection method. Other recommended methods are bw.CvL(), bw.scott() or bw.ppl(). The smoothing kernel used is gaussian, which is the default. Other smoothing methods are: “epanechnikov”, “quartic” or “disc”. The intensity estimate is corrected for edge effect bias by using method described by Jones (1993) and Diggle (2010, equation 18.9). The default is FALSE.
kde_SG_diggle <- density(
childcareSG_ppp,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian") The plot() function of Base R is then used to display the kernel density derived.
plot(kde_SG_diggle)
In the code chunk below, summary() of Base R is used to print the summary report of the
summary(kde_SG_diggle)real-valued pixel image
128 x 128 pixel array (ny, nx)
enclosing rectangle: [2667.538, 55941.94] x [21448.47, 50256.33] units
dimensions of each pixel: 416 x 225.0614 units
Image is defined on a subset of the rectangular grid
Subset area = 669941961.12249 square units
Subset area fraction = 0.437
Pixel values (inside window):
range = [-6.584123e-21, 3.063698e-05]
integral = 1927.788
mean = 2.877545e-06
Before we move on to next section, it is good to know that you can retrieve the bandwidth used to compute the kde layer by using the code chunk below.
bw <- bw.diggle(childcareSG_ppp)
bw sigma
295.9712
Rescalling KDE values
In the code chunk below, rescale.ppp() is used to covert the unit of measurement from meter to kilometer.
childcareSG_ppp_km <- rescale.ppp(
childcareSG_ppp, 1000, "km")Now, we can re-run density() using the resale data set and plot the output kde map.
kde_childcareSG_km <- density(childcareSG_ppp_km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")Next, plot() is used to plot the kde object as shown below.
plot(kde_childcareSG_km)
Working with different automatic badwidth methods
Beside bw.diggle(), there are three other spatstat functions can be used to determine the bandwidth, they are: bw.CvL(), bw.scott(), and bw.ppl().
Let us take a look at the bandwidth return by these automatic bandwidth calculation methods by using the code chunk below.
bw.CvL(childcareSG_ppp_km) sigma
4.357209
bw.scott(childcareSG_ppp_km) sigma.x sigma.y
2.159749 1.396455
bw.ppl(childcareSG_ppp_km) sigma
0.378997
bw.diggle(childcareSG_ppp_km) sigma
0.2959712
Baddeley et. (2016) suggested the use of the bw.ppl() algorithm because past experience shown that it tends to produce the more appropriate values when the pattern consists predominantly of tight clusters. But they also insist that if the purpose of once study is to detect a single tight cluster in the midst of random noise then the bw.diggle() method seems to work best.
The code chunk beow will be used to compare the output of using bw.diggle and bw.ppl methods.
kde_childcareSG.ppl <- density(childcareSG_ppp_km,
sigma=bw.ppl,
edge=TRUE,
kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_childcareSG_km, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")
Working with different kernel methods
By default, the kernel method used in density.ppp() is gaussian. But there are three other options, namely: Epanechnikov, Quartic and Dics.
The code chunk below will be used to compute three more kernel density estimations by using these three kernel function.
par(mfrow=c(2,2))
plot(density(childcareSG_ppp_km,
sigma=0.2959712,
edge=TRUE,
kernel="gaussian"),
main="Gaussian")
plot(density(childcareSG_ppp_km,
sigma=0.2959712,
edge=TRUE,
kernel="epanechnikov"),
main="Epanechnikov")
plot(density(childcareSG_ppp_km,
sigma=0.2959712,
edge=TRUE,
kernel="quartic"),
main="Quartic")
plot(density(childcareSG_ppp_km,
sigma=0.2959712,
edge=TRUE,
kernel="disc"),
main="Disc")
Fixed and Adaptive KDE
Computing KDE by using fixed bandwidth
Next, you will compute a KDE layer by defining a bandwidth of 600 meter. Notice that in the code chunk below, the sigma value used is 0.6. This is because the unit of measurement of childcareSG_ppp_km object is in kilometer, hence the 600m is 0.6km.
kde_childcareSG_fb <- density(childcareSG_ppp_km,
sigma=0.6,
edge=TRUE,
kernel="gaussian")
plot(kde_childcareSG_fb)
Computing KDE by using adaptive bandwidth
Fixed bandwidth method is very sensitive to highly skew distribution of spatial point patterns over geographical units for example urban versus rural. One way to overcome this problem is by using adaptive bandwidth instead.
In this section, you will learn how to derive adaptive kernel density estimation by using density.adaptive() of spatstat.
kde_childcareSG_ab <- adaptive.density(
childcareSG_ppp_km,
method="kernel")
plot(kde_childcareSG_ab)